C     Last change:  ERB  10 Jul 2002    3:08 pm
      SUBROUTINE GWF1EVT6AL(ISUM,ISUMI,LCIEVT,LCEVTR,LCEXDP,LCSURF,NCOL,
     1                  NROW,NEVTOP,IN,IOUT,IEVTCB,IFREFM,NPEVT,IEVTPF)
C
C-----VERSION 14DEC2000 GWF1EVT6AL
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR EVAPOTRANSPIRATION
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
C
C1------IDENTIFY PACKAGE.
      IEVTPF=20
      WRITE(IOUT,1)IN
    1 FORMAT(1X,/1X,'EVT6 -- EVAPOTRANSPIRATION PACKAGE, VERSION 6,',
     1     ' 12/14/2000',/,9X,'INPUT READ FROM UNIT ',I4)
C
C2------READ ET OPTION (NEVTOP) AND UNIT OR FLAG FOR CELL-BY-CELL FLOW
C2------TERMS (IEVTCB).
      CALL URDCOM(IN,IOUT,LINE)
      CALL UPARARRAL(IN,IOUT,LINE,NPEVT)
      IF(IFREFM.EQ.0) THEN
         READ(LINE,'(2I10)') NEVTOP,IEVTCB
      ELSE
         LLOC=1
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NEVTOP,R,IOUT,IN)
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IEVTCB,R,IOUT,IN)
      END IF
C
C3------CHECK TO SEE THAT ET OPTION IS LEGAL.
      IF(NEVTOP.GE.1.AND.NEVTOP.LE.3)GO TO 200
C
C3A-----OPTION IS ILLEGAL -- PRINT A MESSAGE & ABORT SIMULATION.
      WRITE(IOUT,8)
    8 FORMAT(1X,'ILLEGAL ET OPTION CODE. SIMULATION ABORTING')
        call stopfile ! emrl
      STOP
C
C4------OPTION IS LEGAL -- PRINT THE OPTION CODE.
  200 IF(NEVTOP.EQ.1) WRITE(IOUT,201)
  201 FORMAT(1X,'OPTION 1 -- EVAPOTRANSPIRATION FROM TOP LAYER')
      IF(NEVTOP.EQ.2) WRITE(IOUT,202)
  202 FORMAT(1X,'OPTION 2 -- EVAPOTRANSPIRATION FROM ONE SPECIFIED',
     1   ' NODE IN EACH VERTICAL COLUMN')
      IF(NEVTOP.EQ.3) WRITE(IOUT,203)
  203 FORMAT(1X,'OPTION 3 -- EVAPOTRANSPIRATION FROM HIGHEST ACTIVE',
     1   ' NODE IN EACH VERTICAL COLUMN')
C
C5------IF CELL-BY-CELL FLOWS ARE TO BE SAVED, THEN PRINT UNIT NUMBER.
      IF(IEVTCB.GT.0) WRITE(IOUT,204) IEVTCB
  204 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
C
C6------ALLOCATE SPACE FOR THE ARRAYS EVTR, EXDP AND SURF.
      IRK=ISUM
      LCEVTR=ISUM
      ISUM=ISUM+NCOL*NROW
      LCEXDP=ISUM
      ISUM=ISUM+NCOL*NROW
      LCSURF=ISUM
      ISUM=ISUM+NCOL*NROW
C
C7------ALLOCATE SPACE FOR INDICATOR ARRAY(IEVT) EVEN IF ET OPTION IS
C7------NOT 2 or 3, TO AVOID ERROR OF ARRAY (IR) NOT LARGE ENOUGH
      LCIEVT=ISUMI
      ISUMI=ISUMI+NCOL*NROW
C
C8------CALCULATE & PRINT AMOUNT OF SPACE USED BY ET PACKAGE.
      IRK=ISUM-IRK
      WRITE(IOUT,4)IRK
    4 FORMAT(1X,I10,' ELEMENTS IN RX ARRAY ARE USED BY EVT')
      IRK=NCOL*NROW
      WRITE(IOUT,5)IRK
    5 FORMAT(1X,I10,' ELEMENTS IN IR ARRAY ARE USED BY EVT')
C
C9------RETURN.
      RETURN
      END
      SUBROUTINE GWF1EVT6RQ(IN,IOUT,NPEVT,ITERP,INAMLOC)
C
C-----VERSION 20NOV2001 GWF1EVT6RQ
C     ******************************************************************
C     READ EVAPOTRANSPIRATION PARAMETER DEFINITIONS
C     ******************************************************************
C     Modified 11/20/2001 to support parameter instances - ERB
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*4 PTYP
C     ------------------------------------------------------------------
C
C-------READ NAMED PARAMETERS
      IF (ITERP.EQ.1) WRITE(IOUT,5) NPEVT
    5 FORMAT(1X,//1X,I5,' Evapotranspiration parameters')
      IF(NPEVT.GT.0) THEN
         DO 20 K=1,NPEVT
         CALL UPARARRRP(IN,IOUT,N,0,PTYP,ITERP,1,INAMLOC)
         IF(PTYP.NE.'EVT') THEN
            WRITE(IOUT,7)
    7       FORMAT(1X,'Parameter type must be EVT')
        call stopfile ! emrl
            STOP
         END IF
   20    CONTINUE
      END IF
C
C8------RETURN
   60 RETURN
      END
      SUBROUTINE GWF1EVT6RP(NEVTOP,IEVT,EVTR,EXDP,SURF,DELR,DELC,NCOL,
     1     NROW,IN,IOUT,IFREFM,NPEVT,RMLT,IZON,NMLTAR,NZONAR,IEVTPF)
C
C     VERSION 11JAN2000 GWF1EVT6RP
C     ******************************************************************
C     READ EVAPOTRANSPIRATION DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*24 ANAME(4)
      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),
     1          SURF(NCOL,NROW),DELR(NCOL),DELC(NROW),
     2          RMLT(NCOL,NROW,NMLTAR),IZON(NCOL,NROW,NZONAR)
C
      DATA ANAME(1) /'          ET LAYER INDEX'/
      DATA ANAME(2) /'              ET SURFACE'/
      DATA ANAME(3) /' EVAPOTRANSPIRATION RATE'/
      DATA ANAME(4) /'        EXTINCTION DEPTH'/
C     ------------------------------------------------------------------
C
C1------READ FLAGS SHOWING WHETHER DATA IS TO BE REUSED.
      IF(NEVTOP.EQ.2) THEN
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(4I10)') INSURF,INEVTR,INEXDP,INIEVT
         ELSE
            READ(IN,*) INSURF,INEVTR,INEXDP,INIEVT
         END IF
      ELSE
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(3I10)') INSURF,INEVTR,INEXDP
         ELSE
            READ(IN,*) INSURF,INEVTR,INEXDP
         END IF
      END IF
C
C2------TEST INSURF TO SEE WHERE SURFACE ELEVATION COMES FROM.
      IF(INSURF.GE.0)GO TO 32
C
C2A------IF INSURF<0 THEN REUSE SURFACE ARRAY FROM LAST STRESS PERIOD
      WRITE(IOUT,3)
    3 FORMAT(1X,/1X,'REUSING SURF FROM LAST STRESS PERIOD')
      GO TO 35
C
C3-------IF INSURF=>0 THEN CALL MODULE U2DREL TO READ SURFACE.
   32 CALL U2DREL(SURF,ANAME(2),NROW,NCOL,0,IN,IOUT)
C
C4------TEST INEVTR TO SEE WHERE MAX ET RATE COMES FROM.
   35 IF(INEVTR.GE.0)GO TO 37
C
C4A-----IF INEVTR<0 THEN REUSE MAX ET RATE.
      WRITE(IOUT,4)
    4 FORMAT(1X,/1X,'REUSING EVTR FROM LAST STRESS PERIOD')
      GO TO 45
C
C5------IF INEVTR=>0 CALL MODULE U2DREL TO READ MAX ET RATE.
   37 IF(NPEVT.EQ.0) THEN
         CALL U2DREL(EVTR,ANAME(3),NROW,NCOL,0,IN,IOUT)
      ELSE
C    INEVTR is the number of parameters to use this stress period
         CALL PRESET('EVT')
         WRITE(IOUT,33)
   33    FORMAT(1X,///1X,
     1      'EVTR array defined by the following parameters:')
         IF (INEVTR.EQ.0) THEN
           WRITE(IOUT,34)
   34      FORMAT(' ERROR: When parameters are defined for the EVT',
     &     ' Package, at least one parameter',/,' must be specified',
     &     ' each stress period -- STOP EXECUTION (GWF1EVT6RP)')
        call stopfile ! emrl
           STOP
         END IF
         CALL UPARARRSUB2(EVTR,NCOL,NROW,0,INEVTR,IN,IOUT,'EVT',
     1        ANAME(3),'EVT',IEVTPF,RMLT,IZON,NMLTAR,NZONAR)
      END IF
C
C6------MULTIPLY MAX ET RATE BY CELL AREA TO GET VOLUMETRIC RATE
      DO 40 IR=1,NROW
      DO 40 IC=1,NCOL
      EVTR(IC,IR)=EVTR(IC,IR)*DELR(IC)*DELC(IR)
   40 CONTINUE
C
C7------TEST INEXDP TO SEE WHERE EXTINCTION DEPTH COMES FROM
   45 IF(INEXDP.GE.0)GO TO 47
C
C7A------IF INEXDP<0 REUSE EXTINCTION DEPTH FROM LAST STRESS PERIOD
      WRITE(IOUT,5)
    5 FORMAT(1X,/1X,'REUSING EXDP FROM LAST STRESS PERIOD')
      GO TO 48
C
C8-------IF INEXDP=>0 CALL MODULE U2DREL TO READ EXTINCTION DEPTH
   47 CALL U2DREL(EXDP,ANAME(4),NROW,NCOL,0,IN,IOUT)
C
C9------IF OPTION(NEVTOP) IS 2 THEN WE NEED AN INDICATOR ARRAY.
  48  IF(NEVTOP.NE.2)GO TO 50
C
C10------IF INIEVT<0 THEN REUSE LAYER INDICATOR ARRAY.
      IF(INIEVT.GE.0)GO TO 49
      WRITE(IOUT,2)
    2 FORMAT(1X,/1X,'REUSING IEVT FROM LAST STRESS PERIOD')
      GO TO 50
C
C11------IF INIEVT=>0 THEN CALL MODULE U2DINT TO READ INDICATOR ARRAY.
   49 CALL U2DINT(IEVT,ANAME(1),NROW,NCOL,0,IN,IOUT)
C
C12-----RETURN
   50 RETURN
      END
      SUBROUTINE GWF1EVT6FM(NEVTOP,IEVT,EVTR,EXDP,SURF,RHS,HCOF,
     1                  IBOUND,HNEW,NCOL,NROW,NLAY)
C
C-----VERSION 14DEC2000 GWF1EVT6FM
C     ******************************************************************
C        ADD EVAPOTRANSPIRATION TO RHS AND HCOF
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW,HH,SS,XX,DD
      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),
     1          SURF(NCOL,NROW),RHS(NCOL,NROW,NLAY),
     2          HCOF(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     3          HNEW(NCOL,NROW,NLAY)
C     ------------------------------------------------------------------
C
C1------PROCESS EACH HORIZONTAL CELL LOCATION
      DO 10 IR=1,NROW
      DO 10 IC=1,NCOL
C
C2------SET THE LAYER INDEX -- FOR OPTION 1, THE LAYER IS 1;
C2------FOR OPTION 2, THE LAYER IS SPECIFIED IN IEVT.
      IF(NEVTOP.EQ.1) THEN
         IL=1
      ELSE IF(NEVTOP.EQ.2) THEN
         IL=IEVT(IC,IR)
      ELSE
C
C3------FOR OPTION 3, FIND UPPERMOST ACTIVE CELL.
         DO 3 IL=1,NLAY
         IF(IBOUND(IC,IR,IL).NE.0) GO TO 4
    3    CONTINUE
         IL=1
      END IF
C
C4------IF THE CELL IS EXTERNAL IGNORE IT.
    4 IF(IBOUND(IC,IR,IL).LE.0)GO TO 10
      C=EVTR(IC,IR)
      S=SURF(IC,IR)
      SS=S
      HH=HNEW(IC,IR,IL)
C
C5------IF AQUIFER HEAD IS GREATER THAN OR EQUAL TO SURF, ET IS CONSTANT
      IF(HH.LT.SS) GO TO 5
C
C5A-----SUBTRACT -EVTR FROM RHS
      RHS(IC,IR,IL)=RHS(IC,IR,IL) + C
      GO TO 10
C
C6------IF DEPTH TO WATER>=EXTINCTION DEPTH THEN ET IS 0
    5 DD=SS-HH
      X=EXDP(IC,IR)
      XX=X
      IF(DD.GE.XX)GO TO 10
C
C7------LINEAR RANGE. ADD ET TERMS TO BOTH RHS AND HCOF.
      RHS(IC,IR,IL)=RHS(IC,IR,IL)+C-C*S/X
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C/X
   10 CONTINUE
C
C8------RETURN
      RETURN
      END
      SUBROUTINE GWF1EVT6BD(NEVTOP,IEVT,EVTR,EXDP,SURF,IBOUND,HNEW,
     1           NCOL,NROW,NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,
     2           IEVTCB,ICBCFL,BUFF,IOUT,PERTIM,TOTIM)
C-----VERSION 14DEC2000 GWF1EVT6BD
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR EVAPOTRANSPIRATION
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*16 VBNM(MSUM),TEXT
      DOUBLE PRECISION HNEW,RATOUT,QQ,HH,SS,DD,XX,HHCOF,RRHS
      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),
     1          SURF(NCOL,NROW),IBOUND(NCOL,NROW,NLAY),
     2          VBVL(4,MSUM),HNEW(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY)
C
      DATA TEXT /'              ET'/
C     ------------------------------------------------------------------
C
C1------CLEAR THE RATE ACCUMULATOR.
      ZERO=0.
      RATOUT=ZERO
C
C2------SET CELL-BY-CELL BUDGET SAVE FLAG (IBD) AND CLEAR THE BUFFER.
      IBD=0
      IF(IEVTCB.GT.0) IBD=ICBCFL
      DO 2 IL=1,NLAY
      DO 2 IR=1,NROW
      DO 2 IC=1,NCOL
      BUFF(IC,IR,IL)=ZERO
    2 CONTINUE
C
C3------PROCESS EACH HORIZONTAL CELL LOCATION.
      DO 10 IR=1,NROW
      DO 10 IC=1,NCOL
C
C4------SET THE LAYER INDEX -- FOR OPTION 1, THE LAYER IS 1;
C4------FOR OPTION 2, THE LAYER IS SPECIFIED IN IEVT.
      IF(NEVTOP.EQ.1) THEN
         IL=1
      ELSE IF(NEVTOP.EQ.2) THEN
         IL=IEVT(IC,IR)
      ELSE
C
C5------FOR OPTION 3, FIND UPPERMOST ACTIVE CELL.
         DO 3 IL=1,NLAY
         IF(IBOUND(IC,IR,IL).NE.0) GO TO 4
    3    CONTINUE
         IL=1
    4    IEVT(IC,IR)=IL
      END IF
C
C6------IF CELL IS EXTERNAL THEN IGNORE IT.
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 10
      C=EVTR(IC,IR)
      S=SURF(IC,IR)
      SS=S
      HH=HNEW(IC,IR,IL)
C
C7------IF AQUIFER HEAD => SURF,SET Q=MAX ET RATE.
      IF(HH.LT.SS) GO TO 7
      QQ=-C
      GO TO 9
C
C8------IF DEPTH=>EXTINCTION DEPTH, ET IS 0.
    7 X=EXDP(IC,IR)
      XX=X
      DD=SS-HH
      IF(DD.GE.XX)GO TO 10
C
C9------LINEAR RANGE.  Q=-EVTR*(HNEW-(SURF-EXDP))/EXDP, WHICH IS
C9------FORMULATED AS Q= -HNEW*EVTR/EXDP + (EVTR*SURF/EXDP -EVTR).
      HHCOF=-C/X
      RRHS=(C*S/X)-C
      QQ=HH*HHCOF+RRHS
C
C10-----ACCUMULATE TOTAL FLOW RATE.
    9 Q=QQ
      RATOUT=RATOUT-QQ
C
C11-----ADD Q TO BUFFER.
      BUFF(IC,IR,IL)=Q
   10 CONTINUE
C
C12-----IF CELL-BY-CELL FLOW TO BE SAVED, CALL APPROPRIATE UTILITY
C12-----MODULE SAVE THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IEVTCB,BUFF,NCOL,NROW,
     1                          NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV3(KSTP,KPER,TEXT,IEVTCB,BUFF,IEVT,NEVTOP,
     1                   NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
C
C13-----MOVE TOTAL ET RATE INTO VBVL FOR PRINTING BY BAS1OT.
      ROUT=RATOUT
      VBVL(3,MSUM)=ZERO
      VBVL(4,MSUM)=ROUT
C
C14-----ADD ET(ET_RATE TIMES STEP LENGTH) TO VBVL.
      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
C
C15-----MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS1OT.
      VBNM(MSUM)=TEXT
C
C16-----INCREMENT BUDGET TERM COUNTER.
      MSUM=MSUM+1
C
C17-----RETURN.
      RETURN
      END
